#include <fstream>
#include <algorithm>
#include <string>
#include <vector>

#include "sam_parser.hpp"
#include "util.hpp"

SamSequenceMappingParser::SamSequenceMappingParser(const char *fname, std::map<std::string, std::string> &genomes, unsigned int flag) : genome_seqs(genomes)
{
    /*
    char mode[2];
    mode[0] = 'r';
    if(flag & BAM_FILE)
        mode[1] = 'b';
    else
        mode[1] = 0;
    */
    sam_fp = samopen(fname, "r", 0);

    if(sam_fp == NULL)
        error("Failed to open SAM file with 'samopen(" + std::string(fname) + ", \"r\", 0)");

    if(sam_fp->header == NULL)
        error("sam->header is empty!");
    
    /* TODO: figure out how to test for SO:queryname */
    total = 0;
    read = bam_init1();

    status = parse_read();
    last_id = id;
}

bool SamSequenceMappingParser::next_read(SeqRead &sr)
{
    if(!status)
        return false;

    sr.map_locs.clear();
    sr.qual.clear();

    sr.id = id;
    sr.seq = seq;
    sr.qual = qual;
    sr.num = 0;

    if(strand == '-') 
    {
        // strand was negative so the true sequence is not what is given in the
        // SAM file, since that shows it reverse complemented.  In that case
        // reverse it once here before the rest are read.
        std::reverse(sr.seq.begin(), sr.seq.end());
        for(unsigned int k = 0; k < sr.seq.size(); k++)
            sr.seq[k] = rc(sr.seq[k]);
        std::reverse(sr.qual.begin(), sr.qual.end());
    }

    total++;
    
    while(id == last_id && status)
    {
        sr.map_locs.push_back(MapLoci(chr_id, chr_pos, chr_ref, strand));
        status = parse_read();
    }
    
    sr.num = sr.map_locs.size();
    last_id = id;

    /* sr is filled! */
    return true;
}

bool SamSequenceMappingParser::parse_read()
{
    int ret;

    // skip reads that do not map.
    // SAM specification says pos can be 0, but then fails to describe what that means.  Skip them.
    do 
    {
        ret = samread(sam_fp, read);
    } while(ret >= 0 && BAM_FUNMAP & read->core.flag && read->core.pos != 0);

    // TODO: figure out the correct way to test for missing quality scores. 
    /*
    if(bam1_seq(read)[0] == '*')
        error("The file must include the read sequence.");
    */

    if(ret < 0)
        return false;
    
    id      = bam1_qname(read);
    chr_pos = read->core.pos;
    chr_id  = std::string(sam_fp->header->target_name[read->core.tid]);
    strand  = bam1_strand(read) ? '-' : '+';

    qual.resize(read->core.l_qseq);
    seq.resize(read->core.l_qseq);
    
    for(int j = 0; j < read->core.l_qseq; j++)
        seq[j] = bam1_seqi(bam1_seq(read),j);

    for(int j = 0; j < read->core.l_qseq; j++)
    {
        switch(seq[j])
        {
            case 1:  seq[j] = 'A'; break;
            case 2:  seq[j] = 'C'; break;
            case 4:  seq[j] = 'G'; break;
            case 8:  seq[j] = 'T'; break;
            default: error("bad base in sequence.");
        }
    }
   
    // TODO: figure out the correct way to test for missing quality scores. 
    /* 
    if(bam1_qual(read)[0] == '*')
    {
        // do not parse quality, just assign defaults.
        for(int j = 0; j < read->core.l_qseq; j++)
            qual[j] = 1;
    }
    else
    {
    */
        uint8_t *qual_vec = bam1_qual(read);
        for(int j = 0; j < read->core.l_qseq; j++)
            qual[j] = qual_vec[j];
    //}

    chr_ref = genome_seqs[chr_id].substr(chr_pos, seq.size());
        
    if(strand == '-') //reverse
    {
        // in place reverse complement of chr ref.
        std::reverse(chr_ref.begin(), chr_ref.end());
        for(unsigned int k = 0; k < chr_ref.size(); k++)
            chr_ref[k] = rc(chr_ref[k]);
    }

    return true;
}

SamSequenceMappingParser::~SamSequenceMappingParser()
{
    bam_destroy1(read);
    samclose(sam_fp);
}
